Fermion loop simulation of the lattice Gross-Neveu model 
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We present a numerical simulation of the Gross-Neveu model on the lattice using a new rep- 
resentation in terms of fermion loops. In the loop representation all signs due to Pauli statistics 
are eliminated completely and the partition function is a sum over closed loops with only positive 
weights. We demonstrate that the new formulation allows to simulate volumes which are two orders 
of magnitude larger than those accessible with standard methods. 
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I. INTRODUCTION 

Numerical simulations with fermions are notoriously 
difficult. The reason is that the minus signs due to Pauli 
statistics give rise to cancellation effects. In quantum 
field theories the fermions are usually integrated out and 
the fermion determinant appears as a weight factor. Even 
in cases where the fermion determinant is real and posi- 
tive its numerical treatment is very costly, since it essen- 
tially couples all degrees of freedom with each other and 
the individual contributions have changing signs. 

Finding alternative strategies for dealing with fermions 
would considerably improve the quality of numerical sim- 
ulations. Such strategies could either be new algorithms 
(see, e.g., [l| for a prominent example) or a reformulation 
of the problem. Here we discuss the latter: In Q an al- 
ternative representation was given for a two-dimensional 
fermionic quantum field theory, the Gross-Neveu model. 
The partition function was rewritten as a sum over closed 
loops where each contribution has a real positive weight. 
This allows for a new approach to simulate the model 
which avoids dealing with the fermion determinant. 

This paper presents the first test of the loop approach 
for the Gross-Neveu model in a numerical simulation, 
and we explore the prospects and limitations of using 
loop-type representations in a numerical simulation of 
a fermionic system. Our results demonstrate that the 
method is promising and it is worthwhile to pursue it 
in higher dimensions. In higher dimensions loop rep- 
resentations of quantum field theories are known, but 
so far have exclusively been used in the strong coupling 
limit Q ■ Our study here, although 2-dimensional, is per- 
formed at arbitrary coupling. We are currently exploring 
the generalization to higher dimensions and believe that 
for certain four-fermi interactions, representations simi- 
lar to the one used here can be found and successfully 
applied in numerical simulations. 

For gauge theories the situation is more complicated. 
Since gauge fields are oriented one has to use oriented 
loops dressed with the gauge links and complex phases 
appear. This was seen in the Schwinger model, where a 



loop representation exists [J] , but a numerical simulation 
suffers from the fermion sign problem. Upon going to 
the strong coupling limit, the sign problem disappears 
and a numerical simulation with loops again unleashes 
its power [||. Four-fermi interactions may be generated 
with a Hubbard-Stratonovich transformation using scalar 
fields. These do not introduce complex phases and a loop 
representation without signs is possible. 



II. THE LOOP REPRESENTATION 

We begin with discussing the Gross-Neveu model [?J 
and its loop representation. In the continuum the action 
of the Gross-Neveu model is given by S — Sf + Ss with 



Sp — 
Ss = 



d 2 x tp(x) + <p(x) + m] ip(x) , 
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Here ip and ip are Grassmann valued vectors of N fla- 
vors of 2-spinors and we use vector/matrix notation for 
both, spinor and flavor indices. The Euclidean partition 
function is defined by integrating over all fields, 
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Upon integrating out the scalar fields ip, the model turns 
into a purely fermionic theory with a four fermi interac- 
tion given by — g/2 J d 2 x(ip(x)ip{x)) 2 . The Gross-Neveu 
model is well understood analytically (see e.g. 0) and 
has been analyzed on the lattice in various settings @ . 

As it stands, the path integral (J2]) is only formally de- 
fined and a cutoff needs to be introduced. Here we use 
lattice regularization, which replaces the Euclidean space 
time K 2 by a finite regular lattice A. The path integral 
([2]) is well defined when the measure is understood as 
the product over individual measures over the fields liv- 
ing on the lattice points. The action is discretized using 
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the Wilson formulation such that it reads 
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+ f/'C^) + [2 + m]ip(x) 
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For the scalar field ip we use periodic boundary condi- 
tions for both directions, the fermions are periodic in the 
spatial direction and anti-periodic in time. 

Using hopping expansion techniques, the A-flavor lat- 
tice Gross-Neveu model Q can be mapped into a model 
of 2N sets of loops [2]. For convenience we will often 
refer to the loops in different sets as blue, red etc. loops. 
Within each set the loops are non-oriented, closed and 
self- avoiding. However, when loops belong to different 
sets, e.g., a red and a blue loop, they may touch or cross 
each other. The partition function of the lattice Gross- 
Neveu model is then a sum over all possible configura- 
tions of the loops in the 2N sets. Each configuration has 
a positive weight computed from the loops. 

Although [2j gives the partition function for arbitrary 
A, we here only quote the one-flavor expression which we 
use in our simulation. For N = 1 we need two sets of self- 
avoiding loops, red and blue, denoted by r and b. The 
one-flavor partition function in the loop representation 
reads (up to an overall normalization factor) 



z = J2 (iA/2) c(r '% ni(nb V 2 " 2(r 

r,b 



(4) 



In this formula c(r, b) is the total number of corners for 
both red and blue loops. Thus every corner contributes a 
factor of 1/ \[2 to the weight of a configuration. n\ (r, b) is 
the number of lattice sites which are singly occupied by 
either r or b and rii (r, b) is the number of doubly occupied 
sites. We remark that, since the loops in the two sets are 
self avoiding, double occupation can appear only when a 
red and a blue loop cross our run alongside each other. 
The weight factors /i and ji are related to the mass m 
and the coupling g through 

h = (2 + m)[(2 + m) 2 + 5 ]- 1 , f 2 = [(2 + mf + g]' 1 . 

(5) 

We stress that the mapping (H}, is exact in the 
thermodynamic limit. For finite volume different types 
of boundary conditions in the two representations lead 
to finite size effects: In the loop representation we need 
to have closed loops and in a finite volume the loops 
can wind around the compact lattice. The loop configu- 
rations fall into three equivalence classes, C ee , C eo , C°°, 
depending on the numbers of red and blue non-trivially 
winding loops (see also [6(): C ee (even-even): The total 
number of windings for both, red and blue loops is even 
for both directions. C eo (even-odd): One of the colors has 



an odd number of windings for one of the directions. C°° 
(odd-odd) : Both colors have an odd number of windings 
in one of the directions. These equivalence classes cannot 
be linked in a simple way to the boundary conditions in 
the standard representation which we discussed above. 
However, below we will demonstrate that the boundary 
effects vanish as 1/y/V, with V denoting the volume. 



III. NUMERICAL SIMULATION 

The numerical simulation of the Gross-Neveu model 
now is performed directly in the loop representation (|4|) 
using a local Metropolis update (see e.g. (lQ])- We update 
the red and the blue loops alternately, by performing a 
full sweep through the lattice for one of the colors and 
treating the other one as a background field. A sweep 
consists of visiting all plaquettes of the lattice. For each 
plaquette we generate a trial configuration by inverting 
the occupation of the color we want to update for all 
the links in the plaquette. This guarantees that all loops 
remain closed. Furthermore new loops may be generated 
when all links of the plaquette are empty. When the 
trial configuration violates the self-avoiding condition it 
is rejected immediately and the algorithm tries the next 
plaquette. Otherwise the trial configuration is accepted 
with the Metropolis probability p = (l/V2) Ac ff ni f£ n2 , 
where Ac is the change in the number of corners and 
Ani,Ari2 are the changes in the occupation numbers. 
The initial configuration can either be the empty lattice 
(for C ee ) or has one or two winding loops (C e ° and C°°). 

The observables we discuss here are all first and sec- 
ond derivatives of the free energy F — —\nZ, and can 
be written as moments of the occupation numbers. In 
particular for the chiral condensate and its susceptibil- 
ity, which in the standard language are given by 
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we quote the corresponding expressions in terms of occu- 
pation numbers and their fluctuations, 



X 



-p^-(/a(ni) + 2/ 2 (n )) , 

+ [/i-2/ 1 2 / 2 ]((n 1 -(n 1 )) 2 ) 
+ 2/ 2 / 2 ((n + n 1 -(n +?ii)) 2 ) 
-[4/ 1 4 -2/ 2 / 2 ](n )-/ 2 (n 1 ) N ) . (7) 



Here we have introduced no, the total number of empty 
sites, i.e, sites visited by neither a red nor a blue loop. 

Equivalent formulas can be derived for the internal en- 
ergy, the heat capacity as well as for derivatives of the free 
energy with respect to the coupling g. n-point functions 
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FIG. 1: The chiral condensate \ f° r 9 — as a function of 
m for 2 different lattice sizes. We compare the simulation in 
the loop representation (symbols with error bars) to the exact 
result from Fourier transformation (curves). 



may be treated as usually by introducing source fields and 
differentiating with respect to them. This gives expres- 
sions involving correlators of local occupation numbers. 
Finally, the generalization of the above formulas to an 
arbitrary number of flavors is straightforward. 



IV. RESULTS 

In this section we present some selected results which 
serve to illustrate the advantages of the loop approach, 
but also allow to assess its limitations. In order to com- 
pare with traditional methods, we performed a refer- 
ence simulation of the Gross-Neveu model using standard 
methods. The fermions were integrated out giving rise to 
the fermion determinant in a background configuration of 
the scalar field ip. These background configurations were 
computed according to the Gaussian distribution of Ss, 
and the determinant was used as a factor for reweight- 
ing. This is possible, since the eigenvalues of the Dirac 
matrix come in complex conjugate pairs, and the scalar 
field does not have topological modes which could give 
rise to zero eigenvalues. Thus the fermion determinant is 
always strictly positive. 

We stress that the reweighting in the standard formu- 
lation works only in two dimensions due to the numeri- 
cal cost of evaluating the determinant. However, for our 
problem where the scalar fields are independent Gaus- 
sians at each site, reweighting has the big advantage, that 
autocorrelation is avoided. Alternative strategies such as 
Hybrid Monte Carlo, cannot make use of that advantage. 

Another important conceptual point has to be ad- 
dressed: For the free case, g = 0, the standard repre- 
sentation allows for an exact solution with the help of 
Fourier transformation. In the loop formulation, how- 
ever, the case g = is not special at all. Thus g = is the 
optimal point for testing the power of the loop approach 



FIG. 2: Same as Fig. 1, now for the chiral susceptibility C x . 



to the limits because we have exact results on almost ar- 
bitrary large volumes, which we use to compare with the 
data of the loop simulation. Since the weight factors /i 
and f% of Eq. |(5J) are smooth functions of g and m it is 
reasonable to transfer the experience obtained with the 
loop approach at g = to nearby values of g. 

Thus we begin our assessment of the loop approach at 
g = 0. In Figs. 1 and 2 we compare the loop results in the 
C ee sector (symbols) with those from Fourier transforma- 
tion (curves). We use two volumes for the comparison, 
a relatively small lattice of size 32 x 32 and a consid- 
erably larger one, 512 x 512. For the simulation in the 
loop approach at each value of m we typically performed 
10000 sweeps of our local update for both colors for equi- 
libration and used 50000 measurements of the observ- 
ables separated by 10 pairs of sweeps. The observables 
were calculated using the occupation number representa- 
tion (0) and the statistical error was computed with the 
jackknife method. 

Fig. 1 shows that already on the small lattice the data 
points are very close to the exact result. The largest dis- 
crepancy is seen near m = 0, the chiral point where the 
fermions become massless. For the larger lattice the data 
points fall exactly on top of the analytic result. The sit- 
uation is similar for the susceptibility in Fig. 2. For the 
small lattice we find a clear finite size effect, a shift of the 
susceptibility curve. On the larger lattice the agreement 
is almost perfect and only at the chiral point m = we 
still see a slight discrepancy. We remark, that the de- 
crease of the minimum of C x with increasing volume V 
does not signal a phase transition (at N = 1 there is no 
discrete symmetry that could be broken spontaneously). 
The minima decrease only logarithmically with V. For 
g = it can be shown exactly that C x diverges loga- 
rithmically when removing the IR cutoff. For g > we 
could fit the minimum of C x as obtained from the simu- 
lation very reliably with In V. Also the comparison with 
the results from the standard approach shows that the 
largest discrepancy is found near the chiral point, which, 
however, vanishes quickly with increasing volume. 
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FIG. 3: Splitting of the results for the chiral condensate in 
the different equivalence classes as a function of the volume. 

An important part of comparing the standard and the 
loop approach is to test how the different types of bound- 
ary effects scale with the volume and at what rate the 
perfect equivalence of the two representations is reached 
with increasing V. We assess this question directly in the 
loop approach: At a fixed point (m, g) in parameter space 
we compute the chiral condensate x f° r the three differ- 
ent equivalence classes C ee ,C eo ,C°° introduced above. 
This is repeated for several volumes V and in Fig. 3 we 
plot the discrepancy of the results as a function of VV. 
The data shown in the plot are for g = 0.1 and m = 0.0. 
In the plot Ax^ eo ~ ee ^ denotes the discrepancy between 
the C eo and C ee results and Ax (oo_ee) is the splitting 
between C°° and C ee . Also a comparison of \ to the re- 
sults from the standard approach with mixed boundary 
conditions shows a 1/VV behavior. 

The splitting A\ can be analyzed with the help of 
mean field theory. One finds that it should behave as 
Ax oc 1/ y/V. We performed a fit of our data to the form 



Ax oc cV a and found values of a which are close to —1/2 
for all analyzed values of {m,g) (see also Fig. 3). Thus 
mean field arguments as well as our numerical findings 
indicate, that the finite volume effects scale as 1/ y/V. 

V. DISCUSSION 

In this letter we have explored an alternative formu- 
lation for fermionic systems using the example of a 2- 
dimensional quantum field theory. The representation in 
terms of fermion loops allows one to simulate the system 
without having to use fermion determinants. An impor- 
tant aspect is that in the loop formulation used here we 
are not restricted to the case of strong coupling but can 
work at arbitrary g. In this exploratory study we simu- 
late the model using a simple local update and compare 
the outcome to analytic results and the data from a sim- 
ulation in the standard approach. Many observables can 
be expressed in terms of occupation numbers and their 
correlators. We show that finite size effects decrease like 
1 / \/V and thus the thermodynamic limit, where the loop 
representation becomes exact is approached rapidly. 

An important issue is of course the assessment of the 
gain in numerical efficiency when using the loop algo- 
rithm. Already with the local algorithm used here a con- 
siderable increase of the accessible volumes was found. 
Using the same small cluster of PC's the standard ap- 
proach could be used on lattices with a maximum volume 
of 32 x 64, while in the loop formulation we were able to 
simulate systems up to 700 x 700, which is an increase of 
the volume by more than two orders of magnitude. This 
enormous improvement is a strong incentive to search for 
loop representations also in higher dimensional fermion 
systems and for four-fermi interactions no principal ob- 
stacles seem to appear. 
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